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The  governing  differential  equation  for  the 
one-dimensional ,  transonic  flow  in  a  laval  nozzle 
In  the  vicinity  of  the  throat  was  obtained  in  the 
non-diincnsional  form,  A  least  square  finite  ele¬ 
ment  technique  was  used  with  a  linearly  interpol¬ 
ating  polynomial  to  reduce  the  governing  equation 
to  a  system  of  non-linear  algebraic  equations 
which  were  solved  numerically  by  Newton's  method. 
The  system  of  partial  differential  equations  for 
the  two  dimensional  flow  in  a  laval  nozzle  was  also 
obtained  in  the  non-dinensional  form.  The  method 
of  integral  relations  was  used  to  replace  the 
original  system  of  partial  differential  equations 
by  a  system  of  ordinary  differential  equations. 
Using  the  least  square  finite  element  technique  a 
computer  program  was  developed  for  the  construction 
and  solution  of  the  non-linear  equations  for  the 
laval  nozzle  problem.  The  results  including  the 
location  of  the  shock  In  the  flow  are  presented. 


1,  Introduction  \  J 

j 

No  other  aspect  of  transonic  flow  has  been 
the  subject  of  serious  study  for  nearly  so  long  as 
that  of  flows  in  ducts  and  nozzles.  It  Is  remark¬ 
able  that  the  essential  features  of  steady  acceler¬ 
ating  one  dimensional  flows  at  speeds  up  to  and 
beyond  the  speed  of  sound  had  been  well  described 
.  and  understood  by  the  end  of  the  last  century,  be- 
|  fore  the  rapid  development  of  modern  fluid  mechan- 
j ice.  and  a  realistic  analysis  on  a  one  dimensional 
basis  proved  possible.  Xn  particular,  because  v 
;  viscous  effects  In  continuously  accelerating  nozzle 
|  flow’s  are  small,  the  concept  of  boundary  layer  was 
[not  essential  for  progress. 

j  i  ■ 

I  The  survey  papers  of  Hall  and  Sutton  and  of 

!  Sichel ?out line  progress  with  theories  for  nozzle 
flow,  the  latter  paper  in  the  wider  context  of 
viscous  transonic  flow.  Other  newer  techniques  are 
I  due  to  More tti^adap ted  by  Midgal^et  .al.  and 
‘Thompson^.  Holland  others  have  introduced  to  the 
j  nozzle  problem  the  method  of  integral  relations 
which  had  been  successful  in  solving  other  trans¬ 
onic  flow  problems. 

The  problem  of  transonic  flow  In  nozzles  la 
difficult  because  the  governing  equations  are  of 
(he  mixed  type,  changing  from  elliptic  in  the  low 
speed  region  near  the  nozzle  entry  to  hyperbolic  in 
the  supersonic  region.  In  two  dimensional  irrota- 
tional  flow,  the  problem  can  be  solved  by  the  hodo- 
graph  method,  developed  successively  by  Lighthlll?, 
Frankl^  and  Cherry^ In  spite  of  Its  high 
state  of  development,  the  hodograph  method  has 
limitations.  It  can  only  be  applied  to  the  design  ! 
problem  In  plane  flow  and  cannot  be  used  In  the 
direct  problem  of  calculating  the  flow  field  corres¬ 
ponding  to  a  given  contour.  It  does  not  apply  at 


'all  to  axially  symmetric  flow  since  the  presence  of 
'the  radial  term  In  the  equation  of  continuity  des¬ 
troys  the  linear  character  of  the  hodograph  equations. 

For  such  problems  of  more  general  character  it 
is  necessary  to  work  in  the  physical  plane  and  solve 
the  governing  system  of  non-linear  equations  as 
they  stand.  Taylor^used  series  expansion  in  the 
neighborhood  of  the  sonic  throat,  Oswatltsch  and 
Rothstein**developed  an  iterative  technique  and 
Emraons**calculoted  the  floy  in  a  hyperbolic  nozzle 
by  a  relaxation  method. 

A  recent  publication~6nthe  use  of  the  method 
of  integral  relations  for  transonic  flow  in  nozzles 
has  been  published  by  Liddle^1*  ,  Chen*  ,  In  his 
finite  element  solution  to  the  nozzle  problem  was 
unable  to  predict  the  location  of  a  shock  in  the 
flow.  Liddlc 's*^* **method  to  be  valid  requires 
that  there  be  no  shocks  in  the  flow  and  though  this 
method  can  give  solutions  Extending  into  the  super-  J 
sonic  region,  the  solutions  are  invalid  if  there  ara  J 
'shock  waves  or  local  wall  curvature  discontinuities,  q 

!  Tills  paper  sets  out  to  solve  the  direct  problem  I 
of  transonic  flow  In  a  laval  nozzle  and  to  determine  j 
the  location  of  the  weak  shock  using  a  least  square  I 
finite. element  technique^  | 

i 

The  governing  equation  derived  in  the  following  ! 
section  can  be  solved  exactly  by  analytical  methods  « 
and  so  It  was  selected  as  a  model  equation  to  I 

illustrate  the  numerical  method  and  also  because  j 

the  accuracy  of  the  numerical  computation  could  be  i 
compared  with  the  exact  solution,  the  ultimate  aim  .  i 
being  to  use  the  numerical  method  developed  in  this  ( 
paper  for  similar  problems  with  shocks. 


II.  Assumptions  and  Basic  Equations 
For  One-Dimensional  Flow 

Consider  one-dimensional,  isentropic,  steady, 
Inviscld,  transonic  flow  In  a  nozzle  of  varying 
cross-sectional  area  subject  to  no  wall  friction  and 
heat  exchange  with  the  vails.  The  equations  of 
continuity  and  momentum  are 

{  d (pAU)  -  0  (1) 

UdU  +  -  0  (2) 

P 

where,  p  Is  the  density,  A  Is  the  area  of  cross- 
section,  U  is  the  velocity  and  p  is  the  pressure. 

For  lscntroplc  flow  the  speed  of  sound,  a,  is  re¬ 
lated  to  the  compressibility  of  the  fluid  by 
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Eq.  (1)  can  be  written  ns 

UAdp  +  f*!(l'A)  *  0  (5) 

Substituting  for  dp  frc~>  F.q .  (5)  in  Eq.  (4)  and 
using  the  relation  M  -  V/a  where  M  is  the  nach  nua- 
ber,  gives  on  sinplic ir.t  ion 


Finally,  since  this  paper  will  be  prir.arlly  concern¬ 
ed  with  discontinuous  solutions,  it  is  necessary  to 
apply  the  jump  condition 

IV2 J  -  0  (17) 

which  oust  he  satisfied  by  weak  solutions  of  Fq.(14) 


ri(lT)  _  :-J2_  1  <1A  _  . 

dx  2  A  dx 

n  “  J 


(6) 


The  energy  equation  in  differential  form  is 

d(a2  +  V2)  =  0  (7) 


111,  Stat  ir.cn  t  and  Hny  sir  a  1  Cor.i*  cp 
of  the  lava  1  Nozzle  Frobles. 

The  one  diner.sior.al  nozzle  problem  involves  the 
solution  of  the  non-linear  ordinary  differential 
equation,  Eq.  (14)  with  the  boundary  conditions  at 
the  entry  and  the  exit  being  specified  as 


where  y  is  the  ratio  of  specific  heats.  ^DJfferen- 
tiatiug  each  side  of  the  relation  M*  *  (— }  and  sub¬ 
stituting  for  d(aJ)  from  Eq.  (7)  gives  on  simplifi¬ 
cation 


d(U2) 


a2d (M2) _ 

1  +^M2 


(8) 


Substituting  for  d(U  )  from  Eq.  (8)  into  Eq.  (6)  and 
simplifying  yields 


d(H2)  ,  2(1  4  Sj M  ^  1  dA 

dx  ,  u2  A  dx 


<9) 


Assuming  transonic  flow,  let 

M2  -  1  +  6 


(10) 


where  6  «  1.  Assuming  a  parabolic  wall  shape,  let 

2 


1  +  ox 


(ID 


where  a  is  some  constant  depending  on  the  nozzle 
geometry.  Differentiating  Eqs.  (10)  and  (11)  with 
respect  to  x  and  substituting  in  Eq.  (9)  gives 


.  2(1  4  ^  M2)M2 

-  +  - ~ - - - 


2  -1 

Expanding  (1+cix  )  in  a  Taylor  scries  for  small  x, 
i.e.  In  the  vicinity  of  the  throat  and  substituting 
Eq.  (10)  in  Eq.  (12)  and  neglecting  O(x^)  and  0(6*) 
teres,  Eq.  (12)  reduces  to 


V(-a)  -  V.,  V (a)  *  V 

x  n 

where  Vn  is  greater  than  Vi  and  the  jump  condition 
is  as  indicated  by  Eq.  (17). 


Figure  1.  Laval  nozzle. 

Eq.  (14)  is  a  first  order  differential  equation 
which  has  to  satisfy  two  boundary  conditions  and 
thus  it  is  difficult  to  obtain  the  correct  solution 
with  shock  using  standard  numerical  methods.  Sever¬ 
al  numerical  techniques  were  investigated  and  the 
least  square  finite  element  technique  developed  gave 
the  desired  results. 

Solving  Eq.  (14)  analytically  gives 

+  constant  -  ^(py)  *  v  (18) 

For  a  sonic  throat  V  *  0  at  (  *  0,  thus 


*(Y41)oxU  +  «)  -  o 


(13) 


i! 

4 


*"(l~v> 


(19). 


Eq.  (13)  was  scaled  to  give  the  one-dimensional- 
equation  of  transonic  flow  in  a  laval  nozzle  in  the 
normalized  form  aa 
I 

I  -  CO-*)  -  0  (14)  | 


where 

i 


v  “  "  ru  <m2-1> 


Eq.  (IS)  indicates  that  for  M  >  1,  V  is  nega¬ 
tive  and  for  M  <  1,  V  is  positive.  Using  (19),  the 
two  branches  of  the  laval  nozzle  curve  arc  sketched 
In  figure  2.  If  the  flow  In  the  nozzle  was  from 
left  to  right  four  solutions  are  possible,  depending 
on  the  value  of  the  entry  and  exit  Mach  numbers. 

|  Eq.  (14)  can  have  a  subsonic  solution  in  w^iich 
the  flow  accelerates  from  the  inlet  to  the  throat 
and  decelerates  downstream  of  the  throat  and  exits 
at  subsonic  speed  (3-2,  figure  2).  A  continuously 


(13) 


Figure  2.  Two  branches  of  the  Laval  nozzle  curves 

accelerating  flow  solution  can  also  be  obtained  in 
which  the  flow  enters  with  subsonic  speed,  attains 
a  Hach  nur.ber  of  unity  in  the  throat  section  and 
exits  from  the  nozzle  at  supersonic  speed  (3-4, 
figure  2).  A  solution  that  would  be  obtained  for  a 
supersonic  diffuser  in  which  the  flow  could  enter 
with  a  supersonic  velocity,  decelerate  to  reach  the 
sonic  speed  in  the  throat  and  exit  again  with  sub¬ 
sonic  speed  is  the  curve  1-2  shown  in  figure  2.  If 
Vn  >  vl»  the  flow  will  enter  with  subsonic  speed, 
become  sonic  in  the  throat,  attain  supersonic  speed 
downstream  of  the  throat,  have  a  jump  in  velocity 
through  a  shock  and  exit  with  a  subsonic  velocity 
from  the  nozzle*.  The  two  possible  solutions  of 
interest  in  this  paper  are  shown  in  figures  3  and  4. 


e 


,GF  «  U*  satisfying  Hq.  (17).  The  complete  set  of 
Inval  nozzle  curves  arc  plotted  in  figure  6. 


Figure  5.  Graphical  determination  of  the  shock 
solution. 


Figure  6.  Laval  nozzle  curves  for  different  end 
conditions. 


Figure  3.  The  subsonic  solution 


Figure  4.  The  solution  with  shock 


I 


Specifying  Vfl,  the  exit  condition,  the  value 
bf  the  constant  in  Eq.  (18)  can  be  determined  and 
the  shock  location  in  the  model  example  can  be 
obtained  graphically.  Using  Eq.  (19)  AOB  vat  plotted 
as  shown  in  figure  5.  OC  was  then  plotted  as  a 
birr or  Image  of  OB  about  the  £-axis.  Using  Vn,  the 
jronstAnt  in  Eq*  (18)  was  determined  and  the  curve 
pE!U  vas  plotted  on  the  same  figure*  DEHJ  Inter-  \ 
Meets  the  reflected  curve  OC  at  E,  thus  giving  the 
location  of  the  discontinuity.  The  shock  solution  ( 
pi  represented  by  AOCFED  with  the  Jump  condition  _ » 


IV.  Numerical  Procedures  Investigated 

Several  techniques  were  investigated  to  obtain 
the  solution  with  shock  and  the  modified  least 
square  method  was  found  to  give  the  desired  results. 

A  detailed  discussion  of  the  other  methods  investi¬ 
gated  is  beyond  the  scope  of  this  paper,  but  a  brief 
mention  is  made  with  reference  to  them  here. 

The  fourier  sine  transform  technique  was  used 
to  solve  the  governing  equation  and  it  resulted  in 
a  solution  resembling  the  subsonic  solution.  The 
one  and  two  term  least  square  method  using  fourier 
aeries  and  Imposing  the  additional  conditions  that 
V  -  0  (M-l)  and  dV/d£  <  0  at  the  throat  resulted  in 
a  solution  which  did  resemble  the  shock  solution 
but  vas  not  a  good  reproduction  of  the  exact  solu¬ 
tion.  The  results  are  sketched  in  figure  7.  To 
obtain  on  acceptable  result,  a  larger  number  of 
terms  in  the  fourier  series  would  have  to  be  consi¬ 
dered  and  the  effort  required  to  obtain  the  final 

form  of  the  non-linear  algebraic  equations  would  be - 

prohibitive.  Hence  the  simpler  least  square  finite — — 
element  technique  was  preferred  over  the  other 

methods*  - — — -  — 
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1.  Analytical  Solution  with  Shock. 

2.  Fourier  Sine  Transform  Method. 

3.  2-Term  Least  Square  method* 

4.  1-Term  Least  Square  Method 

Figure  7*  Plot  of  V  versus  £  for  two  other  numer¬ 
ical  techniques  investigated. 

V*  The  Least  Square  Finite  Element  Technique 
The  model  equation  considered  was 

-  S(l-V)  -  o  (14) 

with  the  boundary  conditions  V(-a)  «  and  V(a)  « 
Vn.  The  range  between  £  »  -a  to  £  *  a  was  divided 
Into  a  finite  number  of  intervals  and  for  illustra¬ 
tion  purposes  figure  8  shows  4  elements  under  consi¬ 
deration.  Vj  and  Vs  arc  the  specified  entry  and 
exit  conditions  respectively  and  V2.V3  and  V4  were 
computed  as  follows. 


Figure  S.  The  £  axis  divided  into  4  intervals. 

Taking  two  elements  simultaneously  (figure  9) 
and  assuming  a  linear  variation  of  V,  the  system  of 
non-linear  algebraic  equations  was  obtained 
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Figure  9.  Development  of  the  system  of  equations. 

V  -  a  +  b£  {20) 

|trfters  a  and  b  are  functions  of  V  and  £,  that  is  In 
the  interval  to 


Vi-rV/i 


Wi 


.and  similarly  in  the  interval  £^  to 


V  -V 

1  JH 


Let 


R  -  -  C(l-v) 


(21) 


(22) 


(23) 


where ,  R  is  the  residue.  R  *  0  if  the  solution  of 
the  governing  equation  is  exact.  Substituting  for 
V  from  Eq.  (20)  gives 


R  »  2b(a+b£)-(l-a-b£)£ 


(24) 


The  total  error  in  the  two  elements  being  considered 
simultaneously  is  given  by 


ci  2  ci+i  2 

St  Rd£  +  /  Rd£ 

^i-1  M 


Total  error 


(25) 


The  system  of  non-linear  equations  to  be  solved 
was  obtained  by  minimizing  the  total  error  as 


2  ^  2  d  ^  2 

RdC‘° 


(26) 


So*  for  n  unknowns  in  V#  there  are  n  non-linear 
equations  to  be  solved.  Thus  for  three  unknowns  the 
system  of  equations  is 


/f2R2dC  +  3^-  /!3R2dt 

"v2  ^2 


(27) 


wv  ■  w-2 

VWV  -  3^  +  W -  ^3R2dC  ■  °  (28) 

VW  '  3^  'If*  +  3^  'l'*2*  -  0  <29> 

Newton* s  method  was  used  to  solve  the  system 
of  algebraic  equations  and  the  iterative  technique 
is  illustrated  below.  The  (k+l)1^  iterate  is  given 
in  terms  of  the  k-1*1  iterate  as 


jvjfr+D  .  {v}<k>-  [jj"1{F) 


(30) 


where 


»  IJ] 


(J)  is  a  tridiagonal  matrix* 


’.h’A  „ 

”,  ”, 

3F2  3Fj  3F2 

3V^  3V^  3V^ 

0  3V33V4J 

Let 


Ol) 


(X) 


i  .  £ J i { x)  -  -{K}  (-•:> 

s  (X)  vms  obtained  from  Kq.  (32)  by  Caur.fi  ;«n 
tc*l  iclnat  ion  and  vms  then  obtained  from 

Eq.  (31).  The*  initial  values  of  V  chosen  to  Murt 
the  iterative  process  were  arbitrary.  The  conver¬ 
gence  criterion  was 

/»  /  '  v«  l2 
cVi-i  1“tW;  (33) 


for  a  system  with  n  unknowns,  where  C  £  10 

The  solution  obtained  by  the  conventional  least 
square  finite  element  technique  resulted  in  a  solu¬ 
tion  reserbling  the  subsonic  one  shown  in  figure  3, 
The  error  in  each  element  was  computed  by 


f  .  •  ;  .  u  .  of  (*<!. ml  Jems  v.j 
,  .  .  s  was  repeated  for  all  the  c  t  : 

*  rh  set  of  results  thus  obtained  vm  ’ 

tc  the  location  of  thc,.*Jtock.  n\.  i 
step  was  to  check  if  the  Jump  ronditi  :  ,  I 
was  satisfied.  For  example,  cons iih  i  ; :  p  t\i 


-1 

(JO; 

•  hock 


in  e lenient  III,  the  values  obtained  fi  r  V  ,  V  and 

ll  _ _ -I  _ _ r  a _  1  1  ‘  ^ 


were  as  sliown  in  figure  11 


V  -V^i 

V  r  i 


‘v 

- e*1ropoloted  $2 

port  of  the 
curve 


for  each  iteration  and  it  was  observed  that  when 
starting  the  iterative  process  by  using  the  analy¬ 
tically  obtained  values  of  the  unknowns,  the  value 
of  e  was  small  in  all  elements  except  in  the  element 
in  which  the  shock  was  located,  i.e.  e  was  very 
large  in  that  element  (element  IV,  figure  10).  Thus 
the  total  error  would  be  very  large.  The  total  error 


3  win  !  e. 

?r 


Figure  11:  Illustration  of  the  shock  fitting  appro¬ 
ach. 

The  curve  and  V^VS^  arc  the  extrapolated 

parts  of  the  solution.  If  Vgj  -  -  then  it  was 

concluded  that  a  shock  would  occur  in  element  111, 
otherwise  a  shock  v;as  not  possible  in  that  chnent. 
Solution  of  Eqs.  (27),  (28)  and  (29)  by  the  conven¬ 
tional  least  square  method  gave  one  possible  result, 
resembling  the  subsonic  solution,  the  solution  of 
the  same  set  of  equations  by  the  modified  least 
square  method  using  the  shock  fitting  approach  gave 
two  possible  sets  of  results  with  shock,  of  which 
the  solution  with  a  reverse  shock  was  considered  as 
being  physically  not  acceptable.  The  two  existing 
possibilities  of  the  flow  are  shown  in  figures  12a 
and  12b,  for  the  example  considered. 


Figure  10.  Sketch  of  a  discontinuous  solution. 

was  minimized  at  each  iteration  and  finally  the 
result  converged  to  a  solution  resembling  the  sub¬ 
sonic  solution.  So  a  modified  least  square  tech¬ 
nique  was  developed  in  which  the  error  in  the  shock 
element  was  not  considered.  The  method  thus  became 
one  of  searching  and  subsequently  fitting  the  shock 
in  the  proper  element. 


VI.  The  Modified  Least  Square  Finite  Element 
Technique 

The  same  example  as  in  the  previous  section  has 
been  considered  to  illustrate  the  method.  In  the 
modified  least  square  finite  element  technique  the 
shock  is  assumed  to  lie  in  a  certain  element  and  the 
error  in  that  clement  was  not  considered  in  the  solu¬ 
tion  of  the  set  of  equations. 

1  So,  with  reference  to  figure  6,  the  shock  was 
first  assumed  to  be  in  element  11  and  neglecting 

i  /Sv . 


the  system  of  equations,  Eqs.  (27),  (28)  and  (29), 
was  solved  by  Ncvcon*s  method.  Similarly,  assuming 
ho  shock  to  be  In  element  Ill  and  now  neglecting 

.t,  2 


Figure  12a.  Subsonic  solution. 
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j  Figure  12b.  Shock  solution. 

I  To  decide  on  whether  to  accept  the  subsonic 
solution  or  the  shock  solution  as  the  appropriate 
result,  the  total  error  was  computed  for  both  the 
cases  by 


r  ^  r  i  +  lt.2 
L  *  *■  /«  K  d£ 

i-1  1 


which  for  the  cx.vr.ple  considered  becomes 


E  /r  R  d£ 

no  shock  ns 


^  j  1 

+  /r  Vd£  +  /-VdC 


+  //RZd£ 

*4 

9  ^3  o 

E  -  /.  Vd£  +  /f  RZd£ 

s  ^ 

+  /r  R  d£  +  /.  RZd£ 

^3  S 

C5  2 

+  /rVdt 

*4 


Of  the  two  possible  results,  the  solution  with  the 
least  total  error  was  selected,  that  is  if  E^s  <  Es, 
then  a  subsonic  solution  was  the  correct  result  and 
if  Efjs  >  Es,  then  a  shock  solution  was  the  correct 
one. 


o  o 


$hft* 

0  Solution 
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Figure  13.  Plot  of  total  error  versus  element  num¬ 
ber. 


VII.  Results  for  the  One- Dimensional  Problem 


The  governing  equation  (14)  was  solved  numer¬ 
ically  by  the  technique  developed  for  the  set  of 
end  conditions: 


Ytl-O) 


S'Zp-C  Ci-v>«0 

0  6  V('I  0)*0  55122 


Figure  14.  Plot  of  V  versus  £ 
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V  (“0.0737)  c  0  051222 

V  (0*0737)  *  0  0513 


V{-1.0>  -  0.55122  ,  V(1.0)  -  0.6  (38) 

j  V(-0.0737)  *  0.051222  ,  V(0.0737)  -  0.0513  (39) 

j  Both  the  possibilities,  the  subsonic  solution 
jand  solution  with  shock,  were  investigated  and  the 
jproper  solution  was  selected  basing  the  deduction 
!on  the  total  error  in  each  case.  Figures  14  and 

15  are  the  plots  of  V  versus  £  for  the  different 
boundary  conditions.  For  the  boundary  condition  In 
Eq.  (39),  the  values  of  the  Mach  number,  M,  versus 
X  were  plotted  for  a  -  0.25,  indicated  in  figures 

16  and  17. 


Figure  15.  Plot  of  V  versus  £. 


'igure  16.  Distribution  of  the  Mach  number  along 
the  axis  of  a  one  dimensional  laval 
ftoticlc,  a  *  0.1. 
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Figure  17.  Distribution  of  the  Mach  number  along 
the  axis  of  a  one  dir.ensior.al  laval 
nozzle,  a  =  0.25. 


VI 1 1.  A  Note  on  the  Results  of  the  One  Dimensional 
Analysis 

The  results  obtained  in  the  previous  section 
show  that  the  numerical  technique  developed  can 
successfully  be  used  for  similar  problems  with  shocks. 
The  analysis  of  many  flow  problems  assuming  the  flow 
to  be  one  dimensional,  is  a  useful  approximation. 

At  first  sight  this  Is  surprising,  since  the  flow 
model  which  forms  the  basis  of  the  one  dimensional 
method,  differs  profoundly  from  the  actual  physical 
situation.  For  example,  the  one  dimensional  method, 
ignores  the  existence  of  a  non-uniform  and  possibly 
changing  velocity  profile  and  it  assumes  happenings 
at  the  walls  of  the  channel,  e.g.  friction,  to  be 
felt  instantaneously  over  the  whole  cross-section. 
However,  the  one  dimensional  method  was  preferred 
for  its  simplicity.  The  numerically  obtained  results 
and  the  analytical  results  shown  in  figure  6,  are  in 
excellent  agreement.  The  following  sections  are 
devoted  to  the  derivation  and  solution  of  the  two 
dimensional  laval  nozzle  problem. 


IX.  Derivation  of  the  Potential  Equation 

The  nozzle  axis  was  selected  as  the  x-axis  and 
the  origin  of  the  coordinate  system  was  placed  in 
the  center  of  the  throat.  With  the  hypothesis  of 
non. vortical  flow  and  pefect  flow  of  a  perfect  gas 
with  constant  specific  heats  Cp,  Cv  and  y  =  Cp/C^, 
the  potential  equation  for  two  dimensional  flow  is 


3u  ,  2  -2.  3v  ,  2  -2.  0  —  3u 

5T  (.  -«  )  +  ^  (.  -V  )  -  2uv  ^ 


In  this  u,  v  are  the  x  and  y  components  of  velo¬ 
city  and  a,  the  local  sonic  velocity,  is  related  to 
the  critical  velocity  a*  through 

||  .*  -  ^  (u2+v2)  (41) 

» 

\ 

Substituting  Eq.  (41)  in  Eq.  (40),  limiting  the  * 
present  investigation  to  the  vicinity  of  the  throat 
and  Introducing  the  dimensionless  velocity  components 

«.  i 

u  •  eA(l+u)  and  y  -  a*y  (42) 

I 

n  which  u  and  v  are  small  quantities,  the  potential 
equation;  Eq.  (40),  becomes 

>s  +  ♦  nl  v*>  .  *  . _2_  _  ifr-n  u 

h'"  V*  '  »y  Vi  ni 

- =J#J=&-L&  g  Ut-)v_._q - 1 - W 


As  x*0,  y*0;  u  .  ro  ami  a]  -  o  on 

the  basis,  of  Kq.  (42)  ’  .  /  *  i  .  :.st  qui-nl  1  y  tliC 

quotient  v/y  approach  •,  f  he  vclmity  rise 

Du/Dx  along  the  axis  c..(  *.t  i.ot  tr  vanish  at  the 
origin.  Then,  if  the  s  -ill  quantities  u  am!  v  arc 
considered  linear  only,  t!.o  foil*  wing  approx ir.ii c 
relation  is  obtained  f rom  Eq.  (43). 


(Y+l)u  v!  -  +  2v  -  0  (44) 

dx  dy  ey 


Since  the  no;:::le  flow  is  s*..  “  i  t  ri  cal  v:  i : : ; 
respect  to  tl.f  x-axis,  eu/'.y  ;*lso  approach..-  ren 
x>0,  yK).  Consequently  the  urn  2\*  2u/ry  r.jy  be 
ignored  by  means  of  which  llq.  (44)  hecorr.cs 


2y  Y+i  3(u  )  n 

3y  2  3x 


Similarly  the  condition  of  irrotat ionality 
becomes 


1*  _  1“  ■=  o  < 

3x  3y  °  ' 

At  the  vails,  in  the  vicinity  of  the  threat 
for  small  u  and  v  it  can  be  shown  that 


where  (  )'  denotes  differentiation  with  respect  to 


ation  of  the  Method  of  Integral  relations 


The  method  of  integral  relations  was  originated 
by  Dorodnitsyn  in  1958  and  was  subsequent lv  applied 
to  several  fluid  dynamics  probiers .  Several 
scientists  including  Holt,  Liddle  and  Archer  have 
applied  the  MIR  to  the  nozzle  problem.  This  rethod 
provides  a  powerful  tool  for  solving  problems  govern¬ 
ed  by  nonlinear  partial  differential  equations  with 
the  aid  of  computers  and  appears  to  bo  veil  suited 
to  the  laval  nozzle  problem.  It  is  applicable  to 
problems  of  elliptic  or  mixed  elliptic-hyperbolic 
type  in  two  independent  variables.  Certain  depen¬ 
dent  variables  are  represented  as  polynorials  or 
fourier  series  in  one  of  the  independent  variables 
and  the  original  system  of  partial  differential 
equations  is  replaced  by  a  system  of  ordinary  differ¬ 
ential  equations  for  the  coefficients.  These  ordin¬ 
ary  differential  equations  are  then  solved  using  the 
appropriate  boundary  conditions. 

For  the  one  strip  case  the  notations  used  ore 
shown  in  figure  18. 


o.yT* 


Figure  18.  One  strip  notations. 
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n  *  yA 


!k!  a 


x.y  to  x,Ti  v ; . c r c 


jThcn,  as  r)  "  0  on  the  x-nxis  nnd  n  15  1  on  the 
channel  wall  Eqs*  (45) , : (46)  and  (47)  become 


('•  ) 


(v«)n«A-0  <«>' 


*  :  *  ;i  4  Vi!? 

*ljV _ L^_ncit  **1  a  * 

The  me  strip  equations  (57; 

.  V»»llc  wall  shaped  nozzle  redi:. 


u_  -  u  +  a(l -ux2 ) 

l  o 


H_(VJ) 

dt 


t(J-V)  -  0 


(60) 


(f'D 


3v  Y’ 


3v 


3x  “  V  n  3r> 


1  3u 
Y  3jj 


(50) 

(51) 


Integrating  Eqs.  (49)  and  (50)  with  respect  to 
H  between  zero  and  a  general  value  at  a  constant 
x-station  gives 

V?feon-VY+i)^dn 


p  3v  d 

f0  ii  dn 


+  4  Y  <T+1)n  dn  -  0(S2) 


(53) 


For  the  one  strip  case  a  linear  approximation 
for  u  and  v  was  assumed.  For  the  two  strip  case  a 
quadratic  interpolation  could  be  done.  Denoting  the 
Values  on  the  axis  by  suffix  0  and  the  values  on  the 
channel  wall  by  suffix  1, 


■  ■  uo +  *(vV 

v  -  v0  +  n(vrv0) 


(54) 

(55) 


But  Vq  =  0  by  symmetry.  Substituting  Eq.  (51)  in 
Eq.  (55)  gives 


where. 


V  -  - 


2u  +  a(3-ax?) 
o 


ft(l-ax2)  2 _ 

3  '  a(Vl) 


F  in 


2cx_+ E 

2cx2+C 


F  in 


48a3 (y+1) 

2^ 


,  E  =  6a  -  f-q 


(62) 


(63) 


(64) 


G  «  6a  +  /^  ,  c  =  -(y+l)a*  (65) 

q  «=  4c(6  +  a2  (y+1) )  -  36a2  (66) 

with  the  jump  condition  for  weak  solutions  given  by 

0  (67) 


(V2] 


As  indicated  by  the  above  analysis,  it  was 
possible  to  reduce  the  one-strip  equations  to  a 
form  similar  to  the  one  dimensional  governing  equa¬ 
tion  and  this  was  subsequently  solved  by  the  least 
square  finite  element  technique. 

Corresponding  to  the  solution  of  Eq.  (61)  with 
the  boundary  conditions  (39)  as  shown  in  figure  15, 
the  curves  for  uQ  and  uj  arc  plotted  for  a  *=  0.1 
and  a  *  0.25  in  figures  19  and  20. 


nvf 


(56) 


Substituting  for  u  and  v  and  placing  I)  B  1  in 
Eqs.  (52)  and  (53)  gives 


*r  ,u;(2W  +  ui(v2ui” 


.  Y+l  Y*  ,  2  .  2,  2Y* 

+  — r—  S-  |U  +  U  U. -2u,  1 - - 

J  M  ©  Oil  I 


-  »0  +  j  IVY"  -  (Y*)2] 


Assuming  a  parabolic  channel  vail 
I  Y  -  1  +  ox2 


•  0 

(57) 

(58) 

l 

I 

i 

(59) ! 


Figure  19.  Distribution  of  u  along  the  x-axis  and 
>  the  channel  wall  of  a  symmetric  two 

r  dimensional  laval  nozzle  with  parabolic 

arc  vail,  a  ■  0.1. 


where  a  is  a  constant  depending  on  nozzle  geometry, 
(57)  and  (58)  were  solved  for  unknowns  u  and  u. . 

The  two  possibilities  considered  while  reduclngAihese 
equations  further  were  cot2  «  1  nnd  when  order  x2 


terns  were  not  negligible. 
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r..I  Lh-.-  .va  !*s  i  ;  y  a  ;ni  V.  k *. 

been  plotted  in  1 1*«’  cor responding  fif.i.n;. 
suits  clearly  show  iti.it  the*  Jur.p*.  in  vt-Irciiy  <■ 
at  the  fc.ine  value  of  x  for  the  out*  strip  r. 
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Figure  20.  Distribution  of  u  along  the  x-axis  and 
the  channel  vail  of  a  symmetric  tvo 
dimensional  laval  nozzle  with  parabolic 
arc  vails,  a  -  0.25. 

XII.  The  One  Strip  Case  Assuming  ax2  «  1 

Assuming  ax*  «  1,  the  one  strip  equations  for. 
a  parabolic  vail  shaped  nozzle  reduce  to 

u,  *=  u  +  a  (68) 

l  o 
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Figure  21.  Distribution  of  u  along  the  x-axis  and 
the  channel  vail  of  a  symmet r i c  tvo 
dimensional  laval  nozzle  with  parabolic- 
arc  vails,  a  =  0.1. 
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vith  the  jump  condition  for  weak  solutions  given  by 


The  relation  between  M  and  is 
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Figure  22.  Distribution  of  the  Mach  number  elong 
the  x-axis  and  the  channel  vail  of  a 
symmetric  tv:o  dinensional  laval  nozzle 
vith  parabolic-arc  vails,  a  =  0.1. 


and  it  can  be  easily  shown  that 


c  M  —  *1 

n*l  X 

<70 

-  M.  -  1 
*o 

(75) 

The  least  square  technique  was  used  to  solve  the 
governing  differential  equation,  Eq.  (69). 

|  Corresponding  to  the  solution  of  Eq.  (69)  with 
the  same  boundary  conditions  as  In  (39),  the  curves 
•for  u0  and  uj  are  plotted  for  a  ■  0.1  and  a  •  0.25 
‘in  figures  21  and  23  and  the  corresponding  curves 
jfor  Hp  and  are  plotted  in  figures  22  and  24. 

1  An  average  value  of  the  x-vclocity  vas  computed 


average 


**o  +  ! 


(76) 


Figure  23.  Distribution  of  u  along  the  x-axis  and 
the  channel  wall  of  a  symmetric  tvo 
dimensional  laval  nozzle  with  parabolic- 
arc  vails,  a  *  0.25. 

For  a  *  0.1  and  0.25  the  results  indicated  in 
figures  21  and  23  show  good  agreement  vith  the  re¬ 
sults  of  figures  19  and  20  in  the  preceding  section. 
Thus  It  can  be  safely  concluded  that  ox*  <<  1  is  a 
good  approximation  and  it  is  not  necessary  to  go 
through  a  more  complicated  and  tedious  analysis  of 
the  previous  section  to  get  good  results.  The  aver¬ 
age  values  of  the  velocity  plotted  in  figures  21  and 
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Figure  24.  Distribution  of  the  Mach  number  along 
the  x-axis  and  the  channel  wall  of  a 
symmetric  two  dimensional  laval  nozzle 
with  parabolic-arc  walls,  a  =  0.25. 

23  are  interpreted  as  the  average  velocity  of  the 
flow  in  the  channel.  Furthermore  the  shock  occur- 
Ing  at  some  average  value  between  the  wall  and  the 
axis  Is  a  normal  shock  and  the  shock  is  located  at 
the  point  where  there  is  a  jump  in  the  average 
velocity,  (u0  +  a/2).  Figure  25  illustrates  the 
importance  of  the  average  velocity 


1.  ■  ..  + 

6  lor(Y-»l)  2b*  1 


"  =  b*/,8b *7— -«r  * 

V  (ar(Y+l)  2b f) 

Eq.  (79)  was  solved  by  t lie  least  square  technique 
developed . 

Corresponding  lo  the  solution  of  Fq .  (79)  with 
the  boundary  conditions  V (-0.3600)  =  0.23745, 
V(0.3660)  =  0.24250,  the  curves  for  uQ,  uj  and 

“average  and  Mo-  M1  and  Average  arc  Plottcd  in 
figures  26  and  27,  where 


U  -  U  +  ->TT 

average  o  4b* 


Wf>  _M<| 
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Figure  25. 


The  normal  shock  and  the  linear  velocity 
distribution  in  a  symmetric  two  dimen¬ 
sional  laval  nozzle. 
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The  first  study  of  the  continuously-accelera¬ 
ting  flow  in  a  laval  nozzle  was  made  by  Keyer^®. 
The  least  square  finite  element  method  can  also  be 
used  to  obtain  the  solution  of  such  problems.  The 
numerical  and  the  analytical  results  are  compared 
in  Appendix  A. 


XIII.  Application  of  the  Least  Square  Finite  Element 
Method  to  Emmons*  Hyperbolic  Channel  Flow 

The  method  of  integral  relations  Is  now  applied 
to  a  channel  flow  first  calculated  by  Emmons*^.  The 
symmetric  hyperbolic  channel  can  be  described  by 


l  +  *1 

r  1  +  b2 


where,  a  -  0.5646425  and  b  *  0.8253356 

|  The  one  strip  equations,  Eqs.  (57)  and  (58), 
[for  the  Emmons  hyperbolic  nozzle  in  the  vicinity  of 
!the  throat  reduce  to  j 


Figure  26.  Distribution  of  u  along  the  x-axis  and 
the  channel  wall  of  a  symmetric  two 
dimensional  laval  nozzle  with  hyperbo¬ 
lic-arc  walls. 


Fifcure  27.  Distribution  of  the  Mach  number  along 
the  x-axis  and  the  channel  wall  of  a 
symmetric  two  dimensional  laval  nozzle 
with  hypcrbolic-arc  walls. 

I  Emroons^had  solved  the  problem  of  a  two  dimen- 

jeional  flow  of  a  frictionless,  adiabatic,  perfect 
!  gns  Inside  a  hyperbolic  nozzle.  F-mmons*  solution 
was  for  curved  shocks  and  since  the  velocity  field 
-  after  such  shocks  is  not  in  general  irrotatlonal, 
he  had  considered  the  rotation  term  in  the  flow 
following  the  shock  wave.  The  numerical  solution 


-  CU-v)  -  0 


where. 


r.iM  ul?o  I  c  applied  to  tl.ii  u:.‘.tc-.» 
flow  problem. 


.of  I/..-:  era'  u^tuv.  '  \  .  .  :  Ilf  l-ost 

'square  finite  clt  :*■•:. t  i.-3uti*n  for  ti  e  i:.e  strip 
'rase  ihov  an  excel  J «.  .’it  ent.  1  :v  or.  & 1  i,o]ution 

^ls  shown  in  figure  28 
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Appendix  /\ 

I f  the  n or: c! i re r. s i c >r. a  1  velocity  port  *„•  r ha t  i or. s  : j 
and  v  are  defined  by 


Figure  28.  Flow  of  a  compressible  fluid  in  a  hyper¬ 
bolic  channel  as  obtained  by  Emmons. 


(£3) 


In  terras  of  these  conditions  the  i rrotat ional i ty  is 


efu  _  9v 
3y  3x 


(84) 


vhich  enables  a  potential  $  to  be  defined  by 


=  14 

Sx  ’ 


3y  * 


Figure  29.  Distribution  of  rotation;  constant 
rotation  lines;  u)  -  rotation 
w*  -  L'D/ac,  dimensionless  rotation, 

(ac  -  stagnation  acoustic  velocity). 

,  The  flow  downstream  of  the  shock  wave  was  not 
in  agreement  with  Emmons1  solution  due  to  the  fol¬ 
lowing  reasons.  While  Errrons  considered  the  com¬ 
plete  differential  equation  of  compressible  flow  to 
obtain  his  results,  the  least  square  method  was  used 
to  solve  the  approximate  differential  equation  vhich 
gives  excellent  results  for  transonic  flows  in  the 
vicinity  of  the  throat  of  a  laval  nozzle.  As  al¬ 
ready  indicated  before,  only  weak  solutions  were 
investigated  in  this  paper,  while  the  results  ob¬ 
tained  by  Emmons  was  for  strong  shocks.  Secondly, 
while  Emmons  had  considered  rotational  flow  down¬ 
stream  of  the  shock  (figure  29),  this  paper  assumed 
the  flow  to  be  irrotational  before  and  after  the 
shock. 


The  equation  of  continuity  is 


C-24>  -  $ 
x  x 


-  d2  ti 


02)y> 

Y+l  vy/vxx 


Y+l  (1  +  Cx^'y<'xy 


■  ,  2  2ft-l)  .  _-3£l, 2  ,2  - 

5r+l  y+l  be  Y+l  b  "yy 


(S3) 


1 8 

Meyer  in  his  invest igat ion  of  the  continuously- 
accelerating  flow  in  a  laval  nozzle  assured  that  the 
velocity  distribution  along  the  axis  increased  lin¬ 
early,  i.e. 

u  «  Kx,  <f  =  i  Kx2  (S6) 

and  by  direct  substitution  in  Eq .  (85)  of  a  double 
power  series  for  $• 

4  -  l  E  *  xV  ,  (67) 

mm 


*  XIV.  Conclusions 


be  obtained  the  coefficients  6^  up  to  and  including 
the  sixth  order  terms  (cr+n  £  6). 


A  least  square  finite  element  technique  was 
developed  to  solve  problems  of  transonic  flow  with 
shocks.  On  an  average  it  took  less  than  five  se¬ 
conds  of  CPU  time  on  the  IBM  370/Systeras  367  compu¬ 
ter  to  obtain  the  result  for  a  problem  with  nine 
elements  using  this  method.  The  solutions  obtained 
for  the  two  illustrative  examples  show  that  this 
method  is  a  powerful  and  inexpensive  tool  to  solve 
similar  problems  with  shocks.  The  equations  have  to 
be  reduced  in  such  a  manner  that  the  jump  conditions 
can  he  extracted  from  them  and  at  times  this  involves 
Juggling  with  algebraic  quantities,  although  once 
this  is  done  the  solution  is  quite  simple.  The  j 
method  can  be  extended  to  the  two  strip  case  and  | 


4>  f  to2  +  ^Y~  K2xy2  +  ^24-^—  K3y4 


+  K3x2y2  +  ... 


(68) 


It  was  shown  later  that  the  exact  solution  of 
the  approximate  differential  Eq.  (89) 


-  24  6  +  — 4 

vxvxx  y+i  yy 


(89) 


!  i 

jis  the  first  three  terms  of  the  one  given  in  Eq.(88)« 
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o£  the  \  tlo city  i*u  r  t  u:  1  c  t  !^ns  u  .\:>J  v,  Thus  the 
rxa: t  solution  for  the  two  iMt .imv! esuil  transonic 
flow  problem  for  a  t  ont  inuous;  ly-accoloi  at  ing  flow 
becomes 


u  -  Kx  +  KZy2  (90) 

v  ■=  (Y-ll)K2xy  +  K2y3  (91) 


The  equation  of  the  line  along  which  the  flow  is 
parallel  to  the  x-axis  is  given  by 

x  =  -  ^+—  K v2  ( c J 


It  was  also  shown  that  the  length  parameter 
is  related  to  the  radius  of  curvature  Ry  of  the 
streamlines  in  the  throat  region  by 


/(Y+1)R1 

where  Rj  Is  the  radius  of  curvature  of  the  wall  pro¬ 
file  at  the  throat  cade  nondi:r,ensional  by  dividing 
by  the  nozzle  half  height  or  the  radius  at  the 
throat  as  the  case  cay  be. 

With  R^  ■  5,  the  velocity  distribution  along 
the  x-axis  and  the  channel  vails  of  a  symmetric  two 
dicensional  nozzle  with  parabolic-arc  vails  has  been 
plotted  in  figure  30. 


- •  Venn's  solution 
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Figure  30. 


Velocity  distributions  along  the  x-axis 
and  the  channel  vails  of  a  symmetric 
two  dimensional  nozzle  with  parabolic- 
arc  vails,  R.  «  5. 


Equation  (90)  shoves  that  u  varies  quadratically 
in  y  while  in  the  one  strip  case  it  was  assumed 
that  u  varies  linearly  in  y.  Thus  Meyer’s  results 
and  the  results  obtained  numerically  by  the  least 
square  finite  element  method  (LSFEM)  do  not  show 
perfect  agreement.  In  the  two  strip  case  where  u 
and  v  can  be  considered  to  vary  quadratically  in 
y,  while  in  the  Meyer’s  solution  and  the  LSFEM 
solution  are  expected  to  be  in  close  agreement  with 
each  other.  , 
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